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Abstract 

A focused theme in systems biology is to uncover design principles of biological networks, that is, how specific network 
structures yield specific systems properties. For this purpose, we have previously developed a reverse engineering 
procedure to identify network topologies with high likelihood in generating desired systems properties. Our method 
searches the continuous parameter space of an assembly of network topologies, without enumerating individual network 
topologies separately as traditionally done in other reverse engineering procedures. Here we tested this CPSS (continuous 
parameter space search) method on a previously studied problem: the resettable bistability of an Rb-E2F gene network in 
regulating the quiescence-to-proliferation transition of mammalian cells. From a simplified Rb-E2F gene network, we 
identified network topologies responsible for generating resettable bistability. The CPSS-identified topologies are consistent 
with those reported in the previous study based on individual topology search (ITS), demonstrating the effectiveness of the 
CPSS approach. Since the CPSS and ITS searches are based on different mathematical formulations and different algorithms, 
the consistency of the results also helps cross-validate both approaches. A unique advantage of the CPSS approach lies in its 
applicability to biological networks with large numbers of nodes. To aid the application of the CPSS approach to the study 
of other biological systems, we have developed a computer package that is available in Information SI. 
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Introduction 

Systems biology studies how a biological system evolves to 
perform certain function(s), i.e., the design principles of the system 
[1-3]. Reverse engineering is a computational approach to deduce 
biological network structures responsible for given properties [4- 
15]; it addresses situations that while certain dynamic properties of 
a biological system have been observed, the underlying mecha- 
nisms are unknown or incomplete. To reverse engineer, one can 
perform a thorough in silico search to enumerate all possible 
network topologies leading to the dynamic feature at question, 
then identify the most plausible network topologies. To help 
illustrate the basic strategy, let's consider a simple three-node 
network. There are 9 possible links in the full graph including self- 
regulatory links. Each link has three possibilities: activation, 
inhibition, or no presence. Consequendy, there are a total of 3' 
possible network topologies. As a common reverse engineering 
procedure, one can model each network topology against a 



collection of random parameter sets, and evaluate the robustness 
of each network topology (i.e., the proportion of the parameter sets 
allowing each topology to produce the desired dynamic feature). 
This procedure, which we name as ITS (individual topology 
search), has been successfully adopted to analyze several important 
biological processes including segment polarity [7], perfect 
adaption [8], and bistability [13]. However, the ITS approach is 
difficult to apply to systems with large numbers of nodes for two 
reasons. First, the number of network topologies grows dramat- 
ically with the number of nodes. For example, the total numbers of 
network topologies for 4-node and 5-node systems are 
3"' = 4.3xl0' and 3^' = 8.5 x lO", respectively (compared to 
3^ = 2.0x10* for a 3-node system). Second, since the number of 
model parameters increases with the number of nodes, the fraction 
of the parameter space leading to the desired feature (the 'good' 
region) decreases exponentially. For example, suppose that the 
parameter-space dimension is A'^, and that for each parameter, half 
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Figure 1. Flowchart describing tiie CPSS metiiod to determine and analyze network topologies underlying a given dynamic 
property. 

doi:l 0.1 371/journal.pone.01 05833.g001 



of the parameter range falls into the good region, the maximum 
fraction of the good region with the parameter space is 112^ . 

Recently, we developed a computational approach to overcome 
the above difficulties [4] . Instead of modehng individual network 
topologies separately, our approach searches the continuous 
parameter space that defines an assembly of individual topologies 
derived from the full network (with the presence/ absence and the 
strength of each network link indicated by corresponding 
parameter values). Applying this new CPSS (continuous parameter 
space search) approach, we notice that the good regions are 
typically sparse and isolated within the parameter space [4]. To 
efficiendy recover these good regions, we employed a two-stage 
sampling procedure. In the first stage, techniques such as 
importance sample [16] or genetic algorithm [17,18] are used to 
locate isolated good regions with a small number of good 
parameter sets ("seeds") for each region. In the second stage, the 
seeds are used to perform random walk (importance sampling) 
constrained within each isolated good region. This two-stage 
procedure is effective to recover the good regions, even when they 
account for a tiny fraction (e.g., as low as 10~') of the whole 
parameter space. The CPSS method can therefore effectively 



search exponentially growing topological and parameter spatxs 
accompanied with larger networks. For example, we have recendy 
applied this CPSS approach to study a network with 10 nodes and 
64 parameters [19]. 

How well do the results obtained from the two different 
approaches, CPSS and ITS, agree with each other? The present 
work aims to address this question. Here we first present a detailed 
and improved CPSS procedure, we then test the CPSS procedure 
on a previously analyzed problem based on ITS - the resettable 
bistability of an Rb-E2F gene network that controls the 
mammalian cell cycle entr)' [13]. We find that the CPSS and 
ITS approaches, based on different mathematical formulations 
and search algorithms, generate consistent search results. This 
consistency demonstrates the overall effectiveness of reverse 
engineering approaches (CPSS and ITS) in uncovering network 
design principles. To aid the application of the CPSS approach, 
we also develop a computer package that can be conveniently used 
to reverse cngiiKX'r other biological networks, especially those with 
large numbers of nodes with which the ITS approach becomes 
inefficient. 
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Methods 

(a) Problem setup and outline of the CPSS procedure 

The steps below describe our CPSS procedure (with Figure 1 
showing the flow chart). 

(i) Set up a network (with a selected number of nodes) 
associated with the concerned dynamical property, based 
on experimental evidence and literature. 

(ii) Construct mathematical equations to describe the interac- 
tions among nodes in the network. We use the Wilson- 
Cowan function [4,19-21] for its mathematical flexibility. 

(iii) Define the search criteria for the dynamic property under 

consideration {e.g., bistability). 

(iv) Determine the dimension of the parameter space (A^ and 
perform a two-stage search to get 'good' parameter sets 
that fulfill the criteria in (iii). 

(v) Determine the optimal number of clusters formed by the 
good parameter sets. 

(vi) For each cluster, construct the mean value (MV) matrix 
and the coefficient variation (CV) matrix (see detailed 
discus.sions later) to determine "mean" network topolog}' 
and backbone motifs responsible for the desired dynamic 
property. 

(vii) Analyze the roles of individual network links in generating 
the desired dynamic property. 

In this work, the dynamic property we study is the resettable 
bistability of an Rb-E2F gene network that controls the 
mammalian cell cycle entry. Yao et al. demonstrated that the 
Rb-E2F gene network functions as a bistable switch (Figure 2a) 
[13,22]. The Rb-E2F bistable switch converts graded and 
transient serum growth signals into an all-or-none E2F activation, 
which controls the quiescence-to-proliferation transition of mam- 
malian cells. This Rb-E2F bistable switch is resettable: that is, the 
activated switch can be fully shut off when the strength of serum 
signals is reduced below a threshold (point a in Figure 2a) [13,22]. 
The Rb-E2F gene network also exhibits other dynamic properties 
such as the biphasic response of E2F to MYC stimulation [23], 
and likely coordinates the implementation of different dynamic 
properties by constraining associated network parameters [24]. 

To identify the network topology responsible for the emergent 
resettable bistability in the Rb-E2F gene network, Yao et al. [13] 
took a reverse-engineering approach. They coarse-grained the Rb- 
E2F network into a 3-node circuit and examined what network 
topologies among possible link combinations lead to robust 
resettable bistability. Here we follow their setup of the coarse- 
grained network, combining all E2F activators and CycE/cdk2 
into one node EE, all Rb family proteins into another node RP, 
with the third node MD representing the linear signaling cascade 
consisting of Ras, Myc and CycD/cdk4, 6. In this 3-node network, 
9 possible links exist (Figure 2b). We investigate two cases of the 
network. First (Case I), to be consistent with the study of Yao et. al. 
[13], which is based on known network links in the literature, we 
do not consider the link from RP to MD (link 2) and self-regulatory 
links of MD and RP (links 1 and 5, respectively; Figure 2b). 
Second (Case II), we examine the fuU graph that contains aiU the 9 
possible links (see Text SI). 

(b) Mathematical description 

We use the following mathematical formalism of ordinary 
differential equations (ODEs) to describe the 3-node Rb-E2F 
network: 



d[MD] 
dt 



d[RP\ 
dt 



d[EE] 



-yMDiGif^MD-, Wmd) - [MD]) (la) 



Where G((7/; Wj)-- 



= yRp{G{aRp;WRp)-[RP]) 



= yEE{G{(JEE\WEE)-[EE]) 



and 



(lb) 



(Ic) 



Wi 



(Id) 



Here, G( ojy Wj) is a generic ^sigmoidaP function with steepness 
(slope at Wj = 0) that increases with aj for the species. For 
present study, both i andj can be any one of the three nodes, MD, 
RP, or EE. We set a range for aj between 1-10. Wj denotes the 
overall influence of the network on node j. The coefficient (0;; 
measures the strength and direction of the influence of the i 
species on the j*. The term cOjs specifies whether a species j is 
affected by the input signal (serum concentration), and assumes a 
value of 1 for node MD and 0 otherwise. Each ffiij; is a real number 
between [—1, 1], with a positive value for activation and negative 
value for inhibition. Thus, the sign pattern (— , 0, +) of the 'weight' 
matrix cOj; defines the topology of the network. The term cOjo 
determines whether thej* node is 'on' or 'ofF when all other input 
signals are at 0. The parameter Jj determines how quickly each 
species approaches its steady state value. The value of Yj does not 
affect the steady-state behavior of the network but controls the 
non-stationary network dynamics. We denote the concentration of 
a species x by [x]. We construct the ODEs for the 3-node network 
in a formulation similar to previous models [4,20,25,26]. We refer 
[10,27-29] for more detailed discussions and applications of the 
formalism. 

For case I, the model described by Eq. la-d contains 15 
parameters: 6 Wji, 3 yj, 3 aj and 3 Wp. In the following studies we 
fix 7 parameters constant: we set y^E = Ymd — Yrp =1-0 since they 
have no effect on the steady state behavior; we set 
c^MDO — c^EEO — ~0.5 and ojfipQ = 0.5 so that [EE]ss~0 in absence 
of input signals (to be consistent with the experimental observa- 
tion); we also set Oee — 5.0 as a moderate value for the sigmoidicity 
of the output response of node EE. Therefore, our search is in an 
8-dimensional parameter space. For Case II, we consider three 
additional parameters (comdrp, oJmdmd, and coj^p/ip) and perform 
the search in an 1 1 -dimensional parameter space. 

(c) Determination of condition for resettable bistability 

For each gi\(:n parameter set, we calculate the steady state 
concentrations of EE node ([EE]^^) under a set of input serum 
concentrations uniformly distributed between 0 and 10. At each 
input signal level, we use two different sets of initial conditions, 
corresponding to the quiescence state (EE°^^: 
[EE] = [MD] =0.001, [RP] = 1.0 and the proliferation state 
(EE°'^: [EE] = [MD] = 1.0, [RP] =0.001). Temporal evolution of 
each node is simulated using the deterministic improved Euler 
method [30], with sufficiently long simulation time to ensure that 
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Serum concentration 



(a) (b) 

Figure 2. Schematic representation of Rb-E2F networlt. (a) Resettable bistability of the Rb-E2F networl<. Points on the x-axis (a and b) define 
the left and right boundaries of the bistable region, (b) A coarse-grained 3-node Rb-E2F network following the setup of Yao ef. al. [13,22], S: serum 
input; MD: linear signaling cascade consisting of Ras, Myc and CycD/cdk4,6; RP: Rb family proteins; EE: E2F activators and CycE/cdk2. 
doi:10.1371/journal.pone.0105833.g002 



the system reaches steady state. To determine whether a network 
topology creates resettable bistability with a given parameter set, 
we follow the same three criteria used in Yao et al. (that yielded 
top candidate topologies consistent with further experimental tests) 
[13]: 

(i) Being a switch: [EE],,.jn,„ -[EE],,.niin>X. Here, [EE]„.„,a,, 
and [EE] ss-min denote the maximum and minimum values of 
the steady-state EE concentration, respectively, simulated 
with the EE°^^ initial condition for the entire serum input 
range; the threshold X is set to 0.1. 

(ii) Being bistable: A[EE],, = [EE], °^-[EE], °''^>0.1 x([EE],,. 
„i„ - [EE]ss.jnijj) for at least two of the input levels tested 
(except for the saturation level [S]n,ax =10, where the system 
should be monostable). Here, [EE],,"~ and [EE],°^^ 
denote the steady-state EE concentrations obtained with 
the EE*^^ and EE*^^"^ initial conditions under a given input 
level, respectively. 

(iii) Being resettable: A [EE]ss-^0, when [S]^0. 

If all three criteria are met, we refer the associated parameter set 
as a 'good' parameter set. We define a binary score function for 
the k* parameter set: Wj; = — 1 if all the criteria i-iii for resettable 
bistability are fulfilled; Wj^ = 0 otiierwise. 

In addition, to address the concern that the identified candidate 
topologies may be sensitive to the chosen threshold values in the 
above criteria i-iii, we also perform searches with another 
constraint following Yao et al. [13] and others [31]. This extra 
independent constraint, R", is defined as R" = [So£r]""''/[So,J""", 
with [S„fl]™'" and [S„,J™'" correspond to the right and left 
boundaries of the bistable region ih and a. Figure 2a), respectively. 
Thus, R'^ is positively correlated with the width of the bistable 
region. The comparison of the identified candidate topologies with 
or without this extra constraint R" helps evaluate the "robustness" 
of our model predictions. We consider a network topology with 
R''>3.0 as a topology that fulfills the criterion iv; correspondingly, 
we use a score function to define good parameter sets. 



W = S.OxWjj. For convenience of discussion, we refer the searches 
with the criteria i-iii as "NormaP', and those with criteria i-iv as 
"Constrained". Also, all following discussions refer to Case I (see 
subsection a and h of Methods section for definition) unless 
otherwise noted. 

(cO Metropolis sampling algorithm 

Our goal is to sample an 8-dimensional parameter space that is 
bounded and continuous. The sampling algorithm needs to search 
the parameter space thoroughly and generate sample parameter 
sets that are statistically unbiased and significant. Our strategy is a 
random walk method based on the Metropolis Algorithm [32] 
using the following scheme: 

(i) Choose an arbitrary initial 8-dimensional parameter set Qq 
and determine its score: Wo =—1.0 for 'good' set and 
W() = 0.0 otherwise. 

(ii) Generate parameter set 6|;+i from 9|; by Qi,+i = 6k+^C; 

( = 0.025 in this work) specifies the maximum displacement 
per step, and ^ is a vector of random numbers with uniform 
distribution between —1.0 and 1.0. 

(iii) Compute Wj^+i. If W^^Wk+i, accept the step from k to k+ 
1. If Wi;<Wjj+i, accept the step from k to k+1 with 
probabilityp. 

(iv) Update Oj;. If k is larger than a maximum step number 
(Nniax), stop. Otherwise, return to step (ii). 

We implement this strategy in two stages. In Stage I, we set 
p = 0.01, so that the random walk has a tendency to stay in 'good' 
regions of parameter space, but it can also jump out of a good 
region and searches randomly until it falls into another good 
region (or back into the previous region). A random walk with up 
to lO' (Nn^aJ steps is performed in Stage I until a good parameter 
set is identified (^A'k= —1 for the 'Normal' parameter search, and 
W = — 3.0 for 'constrained' parameter search). Under either 
situation {Normal or Constrained), we repeat the Stage I from 
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Figure 3. CPSS software implementation. The major functional modules (above the gray box) are application-independent objects and will not 
need to be changed for a new application. The modules within the gray box are application-specific objects; they inherit all properties of their 
"parent" objects, which help minimize the amount of new coding needed. For example, the present study uses a bistabihty evaluator (under 
"network feature evaluator") and a 3-node system (under "ODE solver"). A different feature evaluator (e.g. adaptation) and a network with larger 
number of nodes can be similarly implemented for a new study. 
doi:1 0.1 371/journal.pone.01 05833.g003 
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Figure 4. Mean network topology underlying the resettable bistability. (a) iVIean value matrix consisting all six network links, (b) Topology 
matrix after discretization of Mean value matrix, (c) Mean network topology obtained from the topology matrix. 
doi:1 0.1 371/journal.pone.01 05833.g004 
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different starting points untU 30-50 'good' parameter sets are 
found. 

In Stage II, eacli 'good' parameter set obtained from Stage I 
serves as a "seed" for performing random walk within the 'good' 
parameter region (with p = 0 to constrain the random walk to tlie 
same good region). The purpose of Stage II is to generate a large 
sample of good parameter sets that occupy each of the different 
good regions within the parameter space. We set N,-,-,3,;= lO' for 
each run and the random walks are sampled every 100 steps, 
generating ~ 1 0"^- 1 0'^ 'good' parameter sets from each seed. 

We note that the introduction of the two-stage Metropolis 
sampling algorithm described above strengthens the efficiency of 
the CPSS method. Stage I sampling helps fmd seeds for 'good 
parameter regions' following random walks. After the detection of 
good seeds in Stage I, Stage II sampling only explores the volume 
element of a 'good region' of each seed in the N-dimensional 
parameter space. This two-stage method directs the search to the 
"functional" (or "good") regions, without spending much time on 
regions not giving rise to the desired properties (resettable 
bistabUity in this case). Therefore, the computational efficiency 
of the CPSS method would be significantly improved over the ITS 
method, especially in high-dimensional systems. 

(e) Clustering of the good parameter sets 

The good parameter sets obtained from the two-stage search 
may form either a single cluster or several different clusters. The 
parameter sets within each cluster correspond to network 
topologies sharing certain common features. Similar to our 
previously developed CPSS method [4], we apply the K-means 



clustering algorithm [33,34] (K = 2 to 12) to identify possible 
clusters of good parameter sets. In addition, here we perform an 
extra procedure to determine the optimal number of K. That is, 
we calculate the mean silhouette coefficient values [35,36] for each 
K; the largest mean silhouette coefficient value suggests the 
optimal K number. If it happens that the optimal K equals 2, we 
would manually check whether the parameter sets are indeed 
distributed in two distinct clusters or they correspond to a single 
cluster (see Text SI for details). 

[f] Discretization of continuous parameter matrix into 
topology matrix 

In order to identify the topological feature underlying the 
resettable bistable mechanism, we coarse grain the continuous 
interaction coefficient ojji (Eq. Id) into a discretized topological 
matrix Tj;. In the topological space, Tj; is only described by (— , 0, +), 
representing inhibition, no interaction, and activation, respectively. 
A cut-off value (0.1) is used to perform the discretization, following 
the rules below: 

(i) If — 0.1^cciji<0.1, Tji is con.sidered as 0 (no interaction). 

(ii) If (Dj;>0.1, Tji is considered as '+' (activation). 

(iii) If (Bj;< — 0.1, Tj; is Considered as ' — ' (inhibition). 



ig) Statistical method to identify the mean motif and a 
backbone motif 

For each cluster identified in step e, we perform the following 
procedure to get a mean motif: 



0.47 

3.12 — 2.63 

0.34 0.56 0.35 

(a) 




Figure 5. Backbone motif underlying resettable bistability. (a) Coefficient Variation matrix, (b) The backbone motif obtained from the 

Coefficient Variation matrix. 

doi:1 0.1 371/journal.pone.01 05833.g005 
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Table 1. Robust minimum motifs underlying resettable bistability for Normal situation. 





Nature of the Link 


Occurrence Probability 


Present as Minimal model In [13]? 


7 9 


98.9 


Yes (as 7 — 9a) 


3 - — 7 


98.1 


Yes (as 2—7) 


3—4—8 


66.0 


Yes (as 2—5—6) 


4—8—9 


66.6 


Yes (as 3—6—9) 


4—6—8 


55.4 


Yes (as 3—5—6) 


6—7—8 


68.6 


Yes (as 5—6—7) 


4—6—8—7 


55.4 


Yes (as 3—5—6—7) 


4—9—8—7 


66.6 


Yes (as 3—6 — 7—9) 


4—8—3—7 


65.9 


Yes (as 2—3—6—7) 


doi:l 0.1 371/journal.pone.Ol 05833.t001 



(i) Calculate the mean of each interaction coefficient (Dj; among 
all good parameter sets. 

(ii) Map the mean value (MV) matrix into a topological matrix 
Tj; using parameter discretization (step f). 

We also define a backbone motif as the motif with the fewest 
number of non-zero 

Oiji that is shared by the maximum number of resettable bistable 
network topologies within a single cluster, and that by itself is 
sufficient to generate resettable bistability. Identification of 
backbone motifs helps to define the core mechanism of resettable 
bistability. To identify a backbone motif, we first calculate the CV 
matrix for each cluster. As the CV matrix measures the dispersion 
of the data along the range of each parameter co^„ a large CV 



value of a given 0)ji suggests that the dispersion around the mean 
interaction strength of the corresponding link is large and thus the 
said link is not essential to be part of a backbone motif Only links 
with CV < Cut-off should be part of a backbone motif For CV > 
Cut-off, we set Tj; = 0 in the backbone motif. To determine an 
optimal value of Cut-off we follow a strategy described as follows. 
As Cut-off decreases, the corresponding motif becomes simpler 
and therefore more network topologies contain this motif; yet, this 
motif cannot be too simple to lose resettable bistability. Therefore, 
there exists an optimal Cut-off value so that the corresponding 
motif has the minimal topology that is sufficient to generate 
resettable bistability, and that the fraction of network topologies 
containing this backbone motif within a cluster is the highest. 
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Figure 6. Pairwise correlation between links 3 and 9: 2D correlation heat map. The x axis denotes tfie link from EE to MD (link 3); the y axis 
denotes the link from EE to itself (link 9). The value on eacfi axis denotes tfie link strength, with tfie positive and negative segments indicating 
activation and repression links, respectively. Color bar on tfie rigfit: tfie fraction of 'good' parameter sets (supporting resettable bistability). 
doi:1 0.1 371/journal.pone.01 05833.g006 
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Figure 7. Pairwise correlation between links 3 and 9: Diagrams of link combinations that correspond to the heat map in Figure 6. 

doi:1 0.1 371/journal.pone.01 05833.g007 



{h) Software implementation 

To help automate our previously developed CPSS procedure 
and improve its computational efficiency, as well as to aid the 
application of CPSS to study a broad range of biological networks, 
we construct an integrated software package with the following 
features. (1) High flexibility: by implementing the CPSS procedure 
in several independent and interchangeable functional modules 
(Figure 3), the need for recoding is minimized to perform CPSS in 
applications involving different network topologies, dynamic 
properties, and different mathematical formulations. (2) User 
friendly: the software is implemented with minimal required 
human-in-the-loop interaction. (3) Computational efficiency: we 
use the C++ programming language to implement the search 
algorithm to minimize execution time, which can perform billions 
of CPSS iteration steps within few days to weeks depending on the 
dimensionality of the problem. An assembled package consisting 
the entire CPSS program is included in Information SI. 



Results 

At stage I, we first search the Rb-E2F ODE system (Eq la-d) 
with 35 arbitrary initial parameter sets, and identify 35 'good' 
parameter sets. For each of these 35 'seeds', we then perform the 
stage II search with lO' Metropolis steps. We obtain a total of 
~2.0xI0^ good parameter sets. We perform the K-means 
clustering (see subsection e in Methods section) and find that all 
good parameter sets form a single cluster in the parameter space 
(See section Tl in the Text SI). 

(a) Structure of the mean motifs and backbone motifs 

Figure 4a gives the mean value (MV) weight matrix denoting 
different link strengths. From the MV matrix we construct discrete 
topology matrix (Figure 4b) following the discretization principle 
(subsection / in Methods section). Figure 4c depicts the mean 
network topology underlying the resettable bistability, which 
contains three positive feedback loops: a positive feedback between 
MD and EE (links 3-7), a double negative feedback between RP 
and EE (links 6—8), and an EE Self-activation (link 9). This mean 
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Figure 8. Mean network motif for resettable bistability under the Constrained situation, (a) Mean value matrix consisting all six network 
links, (b) Topology Matrix after discretization of the Mean value matrix, (c) Mean network topology, which is responsible for resettable bistability 
under the Constrained situation. 
doi:1 0.1 371/journal.pone.01 05833.g008 



network topology (Figure 4c) is consistent with the most robust 
topology in generating a Rb-E2F bistable switch as identified in 
the study of Yao et al. [13]. 

To determine backbone motifs underlying the resettable 
bistability, we construct the CV matrix from the good parameter 
sets (Figure 5a) and set a CV cut-off value 0.36 corresponding to a 
high (more than 0.95) sample ratio. A sample ratio of a given link 
is a fraction of good parameter sets which contain the said link in 
network topologies within a cluster. The identified backbone motif 
(Figure 5b) contains links 7 and 9. That is, this 7—9 motif is a 
minimum circuit that is shared by most resettable bistable network 
topologies; this motif is also sufficient to generate resettable 
bistability by itself This backbone motif is one of the minimum 
motifs discovered by Yao et al. using a different approach, ITS 
[13]. In addition, we identify several robust minimum motifs, 
which are defined as topologies with high occurrence probability 
in the entire good parameter sets. We determine the fraction of 
good parameter sets corresponding to individual minimum motifs 
within the cluster. Table 1 gives a few top minimum motifs with 2, 
3, or 4 hnks. These top minimal motifs are also present amongst 
the top minimum topology hst in the study of Yao et al. [13], 
which signifies the consistency of present results. 

(b) Functional role of the link 3 

Positive feedback regulation is required for generating bistability 
[1]. The mean network topology underlying resettable bistability 
identified in this work (Figure 4c) contains three positive feedback 
loops. Link 3, which is part of the double-positive feedback loop of 
links 3-7, is present in the most robust topology for bistability but 



not in the most robust topology for resettable bistability in Yao et. 
al. [13]. To analyze this discrepancy, we examine the correlation 
between link 3 and other hnks in the mean topology. We found 
within certain parameter ranges, the presence of link 3 (and thus 
the positive feedback loop 3-7) facilitates resettable bistability. For 
example, the heat map in Figure 6 shows that the "good" data 
points are in the upper right quadrant, indicating both links 3 and 
9 should be activation links to support resettable bistability 
(Figure 7). Furthermore, there exists a negative correlation 
between the strengths of links 3 and 9 among the good data 
points; meanwhile, an intermediate strength link 3 (>0.2, <0.6) is 
favorable. These results suggest that the strength of link 3, when 
properly constrained, helps maximize the creation of bistability 
while maintaining its resettability. 

(c) "Lumped parameter" effects 

We note that the above mentioned negative correlation between 
links 3 and 9 to create good data points (Figure 6) provides an 
example of lumped parameter effects. That is, sometimes a 
combinatorial effect of two or more parameters, instead of 
individual parameter values, dictates network dynamics. The 
lumped parameter effects can be explained by the nature of the 
modehng equations (Eq. l(a-d)). Eq. Id explains that the activation 
of the j* species is dependent on the overall net input Wj. As Wj 
combines inputs from all three regulating nodes, any change in 
one parameter, say (BjMDj can be compensated by a change in the 
other two parameters (cOjEE or (Bjrjj, or both). Such parameter 
compensation expands the region of parameter space that supports 
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Figure 9. Backbone motif underlying resettable bistabiiity under tlie Constrained situation, (a) Coefficient Variation matrix, (b) Bacl<bone 

motifs obtained from the Coefficient Variation matrix. See text for details. 

doi:10.1371/journal.pone.0105833.g009 



resettable bistabiiity, and thus enhances the robustness of the 
model. 

(d) Consistent results from the "Constrained" and 
"Normal" situations 

We also perform the CPSS search with the bistable region 
constraint R (see subsection c in Methods section) on the 
resettable bistabiiity. The mean network topology identified 
(Figure 8) shows that the additional constraint (R ^3) does not 
change the mean motif structure obtained in the ^NormaV 
situation. This is also consistent with the observation in [13]. A 



comparison of the MV matrix between the 'Normal' and 
'Constrained' situations suggests that most link strengths in the 
mean network topology under the Constrained situation are 
increased (compared to those under the Normal situation). 
Meanwhile, the backbone motif and the robust minimum motifs 
in the Constrained situation are consistent with those found in the 
Normal situation (Figure 9; Table 2). This consistency suggests 
that the CPSS-identified candidate topologies are not sensitive to 
the chosen threshold values in the search criteria (i-iii, subsection c 
in Methods section). 



Table 2. Robust minimum motifs underlying resettable bistabiiity for Constrained situation. 



Nature of the Link 


Occurrence Probability 


Present as Minimal model in [13]? 


7 9 


99.6 


Yes (as 7 — 9a) 


3 - — 7 


99.0 


Yes (as 2—7) 


3-^-8 


87.3 


Yes (as 2—5—6) 


4—8—9 


87.4 


Yes (as 3—6—9) 


4—6—8 


85.5 


Yes (as 3—5—6) 


6—7—8 


89.2 


Yes (as 5—6—7) 


4—6—8—7 


85.5 


Yes (as 3—5—6—7) 


4—9—8—7 


87.4 


Yes (as 3—6 — 7—9) 


4—8—3—7 


87.3 


Yes (as 2—3—6—7) 


doi:l 0.1 371 /journal.pone.Ol 05833.t002 
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(e) CPSS with all possible links in the Rb-E2F network 
(Case II) 

In Text SI, we also present analysis of Case II, which allows the 
presence of link 1, 2 and 5 (Figure 2b). Notably, results from our 
CPSS analysis posteriori reveal that links 1 and 5 (self-interaction 
links of MD and RP nodes, respectively) are not favorable, 
indicating a functional selection pressure against evolving these 
two links in the Rb-E2F network. On the other hand, the 
additional link 2 (inhibition link from RP to MD) is included in the 
identified mean motif, while the backbone motif obtained from 
CV matrix analysis is unaltered (see section T2 in Text SI). Thus, 
an inhibition from RP to MD (forming a double negative feedback 
loop with link 4, Figure 2b), if evolved in the Rb-E2F network, can 
facilitate the resettable bistability. The link 2, however, is not 
essential for the generation of resettable bistability, as it is not 
present in the backbone motif obtained from the CV matrix 
analysis. 

Discussion 

In this work we present a thorough comparison of two different 
reverse engineering approaches, CPSS and ITS. The ITS 
approach used by Yao el al. [13] examines each possible network 
topologies individually, and adopts the HiU-type function form. 
The present one, CPSS, explores the continuous parameter space, 
and adopts the Wilson-Cowan type function form. Our analysis 
shows that the two approaches give consistent results on the mean 
network topology and the backbone motifs underlying resettable 
bistability in the Rb-E2F network. 

There exist a few quantitative discrepancies between results 
obtained from these two methods, such as the relative ranks (based 
on occurrence probabilities) of the minimum motifs. These 
discrepancies can be accounted by mathematical formalisms. 
Yao et al. use different formalisms to describe multiplicative and 
additive relations among links; we employ a single form (VVilson- 
Cowan) but use variations of parameter values to reflect different 
regulation modes. Also, in their ITS study Yao et al. allow the 
coexistence of two links of opposite signs from one node to another 
node (e.g., link 4 and 5 coexisting from EE to RP, Figure 1 of 
[13]). Unlike the ITS method, the Wilson-Cowan formalism used 
in this work does not treat multiple links between two nodes 
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